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We present a systematical method for obtaining analytical long-range embedded-atom potentials based on the 
lattice-inversion method. The potentials converge faster (exponentially) than Sutton and Chen's power-law po- 
tentials (Philos. Mag. Lett. 61, 2480(1990)). An interesting relationship between the embedded-atom method 
and the universal binding energy equation of Rose et al. (Phys. Rev. B 29, 2963 (1984)) is also pointed out. The 
potentials are tested by calculating the elastic constants, phonon dispersions, phase stabilities, surface properties 
and melting temperatures of the fee transition metals. The results are overall in agreement with experimental or 
available ab initio data. 
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I. INTRODUCTION 

There are quite a few problems in atomistic simulation 
for which long-range potentials are needed. An impor- 
tant one is the problem of structural energy difference 
(SED) . Normally the minimum SED is of the magnitude 
of one percent of the cohesive energy or so. Evidently 
if the potential range is only up to the second nearest 
neighbors, then a pair functional model will predict no 
energy difference between the fee and hep structures. We 
have to extend the ranges of potentials to further dis- 
tance. Usually people impose a cut-off on the potentials 
and adjust the model parameters so that correct SED 
can be produced. However, the unphysical cut-off proce- 
dure thus becomes the dominant factor for predicting the 
SED: Suppose we fix the model parameters and change 
the cut-off distance, then it is highly possible to find that 
the SED varies in sign with respect to the cut-off dis- 
tance(see Fig.l in Ref. 1). The safe way to remove this 
drawback is to extend the range of the potentials so that 
the contributions from the furthest atoms become less 
than one percent of the cohesive energy. Fig. ^ illus- 
trates that to get reliable SED between fee and bcc for 
copper the potential range should be extended into the 
big circle (corresponding to some tolerable error bound). 
Only to that region (and beyond) does the universal bind- 
ing energy relation (UBER) of Rose et al. □ decrease to 
the magnitude of the SED between fee and bcc lattices. 
For alloys the problem will be more complex. There are 
some superstructures with very large unit cells. To cal- 
culate the heats of formation for these competing struc- 
tures needs very long-range potentials. For calculating 
the elastic constants, the potential range is also impor- 
tant. For example, the predicted shear modulus C of bcc 
structure is zero if a nearest-neighbor potential is used, 
so a potential range beyond nearest-neighbor distance 
is required for bcc structure. Long-range potentials are 
also needed for phonon calculation. In the case that the 



unit cell is very large, the potential range should be long 
enough so that all the atoms in the unit cell can interact 
with each other hence the force constants linking them 
do not vanish. 
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FIG. 1. Schematical illustration for the need of long-range 
potential. The filled cirlces denote the radial distribution of 
the atom shells of fee copper. The curve is the UBER for 
copper. The big circle shows the region to which the potential 
range should be extended. 

To determine interatomic potentials, one has to as- 
sume some functional forms for them (such as exponen- 
tial and power-law functions), and the potential param- 
eters are fitted from experimental properties. The pa- 
rameters can be exactly solved within Johnson's nearest- 
neighbor model with exponential potentialsEi and Sutton 
and Chen's long-range model with inverse power-law po- 
tentialsQ. Of course, Johnson's model is not applicable 
in many cases because it is too short-range. On the 
other hand, as inverse power-laws, Sutton and Chen's 
potentials converge, however, slower than the exponen- 
tials (this will be explained later) . In molecular dynamics 
simulation, slow convergence of the potentials will result 
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in the increasing of the time needed to make the neighbor 
hst. In a long-range model that does not employ power- 
law functions, the conventional way to obtain the pa- 
rameters is the numerical method of least-square fitting. 
The disadvantage of the numerical fitting method is that 
it is a little arbitary. Different authors may obtain differ- 
ent parameters, or the same author may obtain different 
parameters at different time, while a small variation of 
the parameters may lead to change of the long-range tail 
remarkable for the SED. This makes the conventional fit- 
ting method problematical when the fitted potentials arc 
to be used to calculate something like phase stability and 
stacking fault energy. These properties are quite sensi- 
tive to the long-range tail which is, however, not very 
refined in the conventional fitting method. 

In this paper, we present a systematical method of 
obtaining analytical long-range potentials with satisfac- 
tory convergence based on the lattice-inversion method 
(LIM). The LIM was first used by Carlsson, Gelatt and 
Ehrenreich (CGE) to get parameter-free nairwise poten- 
tial from ab initio total energy calculation.d Recently, the 
inversion formula of CGE was recasted into a concise for- 
mulism by Chen baS|ed-pn the Mobius inversion transform 
in number theory a^a. Nevertheless, the original two- 
body inversion scheme has, of course, some problems be- 
cause of the lack of many-body contribution. Therefore, 
some manv-body inversion schemes based on the A^-body 
potentialEl and angularly-dependent Stilliiiger- Weber po- 
tential were_developed by Xie and ChenO and Bazant 
and Kaxiras Ej, respectively. 

This paper is organized as follows. In Section II, we 
briefly introduce the LIM. In Section III, we present the 
lattice-iaversion model for the embedded-atom method 
(EAM)Ej. In Section IV, we discuss the parametrization 
procedure in detail. In Section V, we present some cal- 
culated results. The paper is concluded in Section VI. 



II. THE LATTICE-INVERSION METHOD 

The LIM can be traced to an early work by CGE.i The 
idea is to invert a function from its lattice sum which is 
sometimes easier to be obtained. For example, in the 
pair potential model (PPM), the cohesive energy can be 
written as the summation of the pair potential over the 
crystal lattice 
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where Wm is the number of atoms on the m-th shell, Pm 
is the ratio of the radius of the m-th shell to the nearest- 
neighbor distance. The cohesive energy as a function of 
lattice spacing can be calculated by using first-principles 
method, or simply taken as the UBER. Then as CGE 
suggested, one can use the following inversion formula to 
obtain the so-called ab initio pair potential V{r) 
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The multiple summations make the ordering for the in- 
version coefficients not obvious. It was not until recently 
that Chen put forward his elegant Mobius inversion for- 
mula on three-dimensional crystals E The Chen-Mobius 
inversion formula is very simple 



V{r) = 2 ^ timE{prar) 



(3) 



The Mobius coefficients /im can be determined by 
^1 l/wi 

fj.,n = -{l/wi) ^ pikWi (m > 2) (4) 

Pk \Pm,k^m 

where I is the natural number which satisfies pi = Pm/Pk- 
Obviously, Chen's formula requires the set P = {pm\'m G 
N} should be a multiplication-close one, i.e., given two 
arbitrary elements Pi,Pj G P, their product piPj should 
be in P too. Actually the crystal lattices sc, bcc, fee, hep 
and diamond etc. do not automatically satisfy this re- 
quirement. Therefore, before applying the Chen-Mobius 
formula we have to at first construct a close set Q, 
which should include at least part of the elements in the 
orginal set P. This task is easily done for sc and fee: 
For them the set P is simply + j"^ + k^Ri \i, j, k e 

Z,i'^ + -f A:^ _^ ^g^^ construct a new set 

Q={v^-Ri|^ e N}, which covers P. The numbers of 
atoms on the shell ^/nRi vanish if n cannot be written 
as the square sum of three natural numbers -I- + . 
But for other lattices such as bcc, it is difficult to find 
a natural close set which covers all the elements in 
the set P= {^Ji^ +p + k'^a\i,j,k G Z,i'^ + f + k'^ ^ 
0} U V(« + 1/2)2 + + 1/2)2 + (fc + l/2)2a|i, J, fc e Z}, 
where a is the lattice constant. However, for the present 
physical problem we do not have to construct a close set 
covering all the elements. Note that the expansion of 
eq.(j^) should be convergent. That is to say, usually we 
can truncate at some shell, say, the M-th shell, beyond 
which the function V{r) has become small enough to be 
neglected. So we can approximate eq.(0) by E{Ri) — 
12m=i (PmRi) , then we have a set with M elements 
P= {pi,p2, ■ ■ ■ ,Pm}- We can easily generate a close set 
Q which covers P: Q= {p^Pa' ' ' 'Pm' l^i, fe. ■ ■ ■ .ku ^ 
0, 1, 2, 3, • • • , fci -f ^2 + • • • + kM ^ 0}. Re-ordering this 
close set from the smallest element to the biggest one, we 
get the set as Q= {Prn\Pm < -Pm+i:"^ G N}. Then we 
can rewrite eq.(|l|) as E{Ri) = (1/2) X^^^i 
where — Wm when p„ — p„i and vanishes when p„ 
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equals none of the elements in P. The inversion is sim- 
ply the same as eq.(§), with the Mobius coefficients /i„ 
determined from w„. 

In the one-dimensional case, eq.(^ becomes the 
number-theoretic Mobius inversion formula 



V{r) 



fi(n)E{nr) 
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where fi{n) is the number-theoretic Mobius function. 
This jjphysical mapping (Fig.||) was first dicovered by 
Chenl. 




FIG. 2. A physical mapping of the Mobius inversion theo- 
rem in arithmetic number theory. 

Very recently, Bazant and Kaxiras have presented 
a novel scheme to obtain effective angularly-dependent 
many-body potentials for covalent materials by inversion 
of cohesive energy curves from many configurationsEj. 
Theiij-ipethod, together with our previous work for the 
EAMIlil, have shown the potentiality of the LIM as a 
shortcut to obtain the effective interatomic potentials. 



III. THE LATTICE-INVERSION 
EMBEDDED-ATOM MODEL 

Generally the physical properties we use to fit the po- 
tential parameters are related more closely to the lat- 
tice sums than to the individual potential function them- 
sleves. Given the lattice sums of the pair potential V{r) 
and electron density /(r) 
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We can easily rewrite the formulas of the bulk modulus 
and Voigt shear modulus 

i 

+ 2F" {pe)[Y,R^f' {R^)f} (») 



(9) 



as 
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where is the atomic volume, T4ff (r) is the effective pair 
potential Ves(r) = V{r) + 2F'(p)/(r). And for trans- 



forming eqs. dq) and ( 
forms of eqs.pOl) and 



to the effective nearest neighbor 
, the following relationship has 
been used: if there is a function h[r) whose lattice sum 
is another function H{Ri), then 



= RiH^^'HRi) (12) 
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Thus we have an equation related to the Cauchy pressure 
9Bn~15Gn = F"{p)Rl[p'iRi)f (13) 
On the other hand, the vacancy-formation energy 



E, = -<P + J2{F{p - fiRi)) - F{p)] + E, 
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can be approximately written as the lattice sum of the 
effective pair potential E^ = (1/2) VcsiRi), since the 
numbers of atoms on the shells are much greater than 
1, and the negative relaxation energy further reduces the 
error. This approximation has been checked by a simple 
relaxation calculation in which only the nearest neighbor 
atoms around the vacancy site are allowed to relax. It is 
found that in the case of copper the calculated unrelaxed 
vacancy-formation energy is 1.34 eV, while the relaxed 
result is 1.31 eV, closer to the experimental value. Hence, 
the difference between the vacancy-formation energy and 
the sublimation energy can be written as 



Es^E.^F {p)p{Ri)^F[p{Ri) 
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We can see from eqs.(y_3|) and (|15|) that the nonlinearity 
of the embedding function reflects the many-body nature 
of the embedded-atom potential. If F{p) is a linear func- 
tion with respect to p (corresponding to a PPM), then 
we have 3B = 5G (the Cauchy relation) and Eg = Ey, 
which are the two well-known drawbacks of the PPM. 

The nearest-neighbor distance of the equilibrium lat- 
tice can be obtained by minimizing the total binding en- 
ergy 
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Another condition we need to consider is the normali- 
sation for the electron density. Integrating both sides of 
eq-(0) with respect to Ri 
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Note that the electron density f{r) should be normalized 
inr^ f {r)dr = N (where N is the number of elec- 
trons), we obtain 



4nRlp{Ri)dRi = S{3)N 



(18) 



where S{3) = J2m^™./Pm- alloy case, the parameter 
N should be determined by considering the charge trans- 
fer. This consideration is based on empirical Miedema's 
equation which well decribes the heats formation of bi- 
nary alloys. The attractive term in Miedema's equation is 
re-interpreted by Pettifor as the contribution of the elec- 
tronegativity difference, which is related to the charge 
transfer Jlil 

The embedding function in the present model is as- 
sumed to be a power-law one 



F{p) = -Ap 



i/x 



(19) 



A = 1 corresponds to the PPM, while A = 2 corijesponds 
to the TV-body potential of Finnis and Sinclair.B In the 
following section we shall show that with appropriate 
functional forms for the electron density and pair po- 
tential, this embedding function will produce exactly the 
UBER, and the parameter A is insensitive to the func- 
tional forms of the electron density and pair potential. 

The electron density and pair potential in the present 
model are structure-dependent as we can see from their 
inverted formulas 

/(r) ^ pip{pir) + P2p{P2r) + PsPiPsr) H (20) 

V{r) ^ 2/ii$(pir) + 2^2$(p2r) + 2fi3<P{p3r) + ■■■ (21) 

The functions of /(r) and V{r) are the linear combina- 
tions of their lattice-summed functions p{R) and $(i?), 
while the structural dependence is included in the Mobius 
inversion coefficients pm and the radius ratio Differ- 
ent kinds of functions for p{R) and <&(-R) will be used to 
control the potential convergence, as shown in the next 
section. 



IV. PARAMETRIZATION 

A. p{R) and $(-R) are exponential functions 

It has been found by Banerjea and Smith using the 
effective-medium theory that the off-site electron density 
exhibits a universal relationship with respect to lattice 
spacing: p* = exp(— a*), which was used to explain the 
physical origin of the UBEjE. within the framework of 
local density approximationllj. Based on the results of 
Hartree-Fock calculations, Mei, Davenport and Fernando 
also pointed out that the lattice sum of the electron den- 
sity as a fuHction of lattice constant shows exponential 
behaviour .E3 Therefore, it is plausible to take p{R) as an 
exponential 



p{Ri) = PeCXp 
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As a short comment, pHie would like to point out that when 
the authors of Ref. t£l came to the above equation they 
just used a complex function as f{r)=fe X]f=o ci{Rie/ry 
to fit it. One can see in the present method we do not 
have to fit. The individual function is accurately given 
byeq.(|§). 

The repulsive energy is often assumed to have a rela- 
tion with the bond energy (i.e. the embedding energy 
in this case) like Uicp{R) oc [Uhond{RW , where 7 is 2 
accordins-io the so-called Wolfsberg-Helmholtz approxi- 
mation. 113 Therefore, we assume $(-R) is also an expo- 
nential function 



<i>(i?i) = <i>e cxp 
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In our method the parameters can be exactly solved as 
if the model were a nearest neighbor one. The solutions 
are 



A = 



/3 = 
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(25) 
(26) 
(27) 
(28) 
(29) 



The binding energy equation is a Morse-like function 
£^coh(i?i) = *eexp 

-Api/^exp -:^(^-l) (30) 
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Different from other equations of state, eq. (pOD includes 
the inputs of the Cauchy pressure and vacancy-formation 
energy. 



B. p{R) and $(i?) are gaussian functions 

It has been shown that in the present method all 
the parameters are analytically determined by the input 
physical properties, which are only for the equilibrium 
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lattice. The potential convergence depends on the func- 
tional forms we take for p{R) and $(-R)- The gaussian 
function is an alternative choice 



= peexp 
$(i?i) = $eexp 
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(31) 



(32) 



The solutions are the same with the above subsection 
except 



(33) 
(34) 




The binding energy equation is 



EcohiRl) = *ecxp 
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Apl^^ exp 
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C. p{R) and ^{R) are modified exponential functions 

The electron density may not be a simple exponential 
function. As we know it is always a combination of some 
Slater orbitals r^e^^^. In order to reflect this, we suppose 



p{Ri) = ( ) exp 
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The pair potential remains the same as eg. (p^) . The 
solutions of A are 
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1 + An + ?1 



Es — Ey 

9Bn 
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where Aq — 5GEs/3BEy. In this case, we have two solu- 
tions. Each of them can exactly reproduce the physical 
inputs. This simple example then implies that there may 
exist several different attractors leading to different re- 
sults when the conventional fitting procedure is used to 
search for an approximate solution. This problem may 
merit a thorough investigation, and will not be discussed 
in the present paper. Note that {Eg — Ey)/18Bil, ^ 1 



and {Es — E.i,)/15G^ ^ 1, if n is taken to be 1, then the 
approximate solutions will be A"*" = Aq and A~ = 1. The 
latter solution is just the PPM which is then excluded. 
The solutions for the remaining parameters are 
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(38) 
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The expressions for the other two parameter <I>e and A 
are identical with those presented in the first subsection. 
The binding energy equation is 



^^coh(-Ri) = ^-eexp 
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When a/X = P and A = n, the above equation is just the 
UBER 
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From the above subsections, one can find two points 
to support the power-law embedding function. The first 
point is that this embedding function (given by Apl^'^ 
and A) is independent on the given funcional forms of 
the electron density and pair potential. This physically 
underpins the local nature of the embedding function. It 
is also true when the pair potential and electron density 
take the power-law forms. The second is that the given 
binding energy equation is very close (and even identical) 
to the UBER. This consistency is necessary for a good 
description for the thermal expansion (the anharmonity 
effect )E3. In some sense, the method may represent an 
embedded-atom explanation for the UBER. 

If p{Ri) and ^(-Ri) take the form of power law the 
solutions for the parameters will remain unchanged ex- 
cept that Pe cannot be determined since the power- 
law function cannot be normalized. However, an al- 
(37) ternative parameter Apl^^ can be determined, as 



is the case of Sutton-Chen's potential. The inverted 
functions, for example the pair potential, is given as 
V{r) = ^e/S{(5){r/Rxeyl^. While in the case of expo- 
nential, V{r) < ($e/12)e'^[exp(r/Eie)]-''. Since S{(}) is 
only a little greater than 12 and the paramters $e, [3 are 
the same in the two cases, the above two functions cross 
approximately at the NN distance. Beyond the NN dis- 
tance, the exponential potential is smaller and decreases 
much faster than the power law. 
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A. Structural stabilities 
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FIG. 3. The inverted electron densities and pair potentials 
(inset) for the noble metals. The unit of /(r) is A~^. 

Fig. 1^ shows the typical shapes of the inverted pair 
potential and electron density, respectively. One can see 
from the figure that the inverted electron densities and 
pair potentials decrease rapidly. In our calculation, the 
cut-off is placed at about 3ae, where the pair potential 
and electron density have been negligible. 

V. APPLICATIONS OF THE POTENTIALS 

The physical inputs for the present model are listed in 
Tab. The elastic moduli of Al are from Simmons and 
Wangij, those of Ni, Pd^-Pt, Cu, Ag and Au are from 
Foiles, Baskes and DawsE^I; Lattice constants and cohe- 
sive energies are all from KittelEd; Vacancy-formation en- 
ergies foE_fcc transition metals are from Eoiles, Baskes 
and DawEl, and that for Al is from BallufiB. We do not 
list the number of electrons since Pe can be incorparated 
with ^ as a parameter it is not used in monoatomic cal- 
culations. It is only important in alloy calculations, in 
which it descibes the charge transfer effect. 



TABLE I 


The model inputs 


aB,Ea,Ev, B,G. Oe is in 


A, Es 


and Ev are 


in eV, B and G are 


in 10" 


N/m^ 




Element 




Es 


Ev 


B 


G 


Ni 


3.52 


4.44 


1.7 


1.804 


0.95 


Pd 


3.89 


3.89 


1.54 


1.95 


0.54 


Pt 


3.92 


5.84 


1.6 


2.83 


0.65 


Cu 


3.61 


3.49 


1.3 


1.38 


0.55 


Ag 


4.09 


2.95 


1.1 


1.04 


0.34 


Au 


4.08 


3.81 


0.9 


1.67 


0.52 


Al 


4.05 


3.39 


0.7 


0.76 


0.266 



Phase stability is the first test for the long-range po- 
tentials. For copper, the EAM result i?fcc-bcc (25.8meV) 
falls in the middle of the nonrelativistic and semirela- 
tivistic ab initio values (-17.7 ipeV and -48.8 meV) re- 
ported by Lu, Wei and ZungerE^. The EAM result for 
i'fcc-diamond cquals 1.07 eV, also close to their ab initio 
result (1.35 eV) for diamond-like copper with the correc- 
tion of nonspherical charge-density inside the muffin-tin 
sphere. For all the studied elements, the present model 
predicts the fee structure to be the ground state (see 
Tab. |lu). However, similiar to Sutton and Chen's po- 
tentialsQ, the EAM result for i?fcc-hcp is virtually zero, 
so we did not print the binding energy curve for hep in 
Fig.^. This failure is believed to be due to the absence 
of the angularly-dependent or higher order moment con- 
tributions. It has been pointed out by Ducastelle and 
Cyrot-Lackmann that it is mainly the third and fourth 
moments that are respeasible for the SEDs among bcc, 
fee and hep candidates.Eil 



v,\ \ 


Cu 


^V, \ diamond 




^ \ \ 




VA _sc_.,---' 




\ \ 




bcc 


fcc^ 







2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 
Nearest neighbor distance (angstrom) 



FIG. 4. The predicted phase stability of copper using the 
present model. 



TABLE 11. The predicted structural energy differences for 
the fee metals. The energies are in eV. The numbers in the 
parentheses denote the function types used: l=exponential; 



2=gaussian 


; 3=modified 


exponential. 




Element 


Etcc — Esc 


Eicc Elccc 


Efcc -^'diamond 


Ni(2) 


-0.58 


-5.31x10"^ 


-1.36 


Pd(l) 


-0.53 


-3.42x10"^ 


-1.21 


Pt(l) 


-0.62 


-4.67x10"^ 


-1.39 


Cu(l) 


-0.44 


-2.58x10"^ 


-1.07 


Ag(l) 


-0.40 


-2.28x10"^ 


-0.94 


Au(l) 


-0.35 


-2.29x10"^ 


-0.83 


Al(3) 


-0.28 


-2.11x10"^ 


-0.66 
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B. Elastic constants and phonon eigenfrequencies 



The comparison of the calculated and experimen- 
tal data for the elastic constants Cn, C12, C44, the 
anisotropy ratios C/C" (C = C44, C = (Cn - Ci2)/2), 
and the phonon longitudinal and transverse frequencies 
7l, 7t at the boundary of the Brillouin zone are shown 
in Fig. ||. The elastic constants were calculated by exert- 
ing the corresponding strain matrices to the lattice. The 
vibrational eigenfrequencies are calculated by diagonal- 
izing the EAM dynamical matrixE^. The experimental 
data for the elastic constants of Ni, Ed, Pt, Cu, Ag and 
Au are from Foiles, Baskes-and Dawo, those for Al are 
from Simmons and WangEj. The experinierital data for 
the phonon eigenfrequencies are from RefBj. The calcu- 
lated results are overall in agreement with the experimen- 
tal data. Fig. ^ shows the predicted phonon dispersion 
curves for copper along the high symmetry directions are 
in excellent agreement with the experimental data. 
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FIG. 6. Comparison of the theoretical phonon dispersion 
curves of Cu (solid lines) with the experimental data (filled 
circles) along the high symmetry directions. 

It has to be pointed out that when the present 
model is applied to the bcc transition metals with low 
anisotropic ratios, the calculated elastic constants C" and 
C44 severely disagree with the experimental data, despite 
of that the bulk modulus and the Voigt shear modulus 
can be reproduced. This may imply that for the bcc 
transition metals the directional bonding is significant. 
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FIG. 5. Comparison of calculated and experimental re- 
sults for the elastic constants, the longitudinal and transverse 
phonon eigenfrequencies at the Brillouin-zone boundary, (a) 
Cn; (b) C12; (c) C44; (d)7i; (e)7T. 



C. Surface properties 

To calculate surface properties, we employ a simulation 
box with size lOx lOx 10 and periodically reproduced in 
X,?/ and half z directions (rather than a slab). For (111) 
surface, a periodic boundary condition with rhombic ge- 
ometry has been applied. Relaxation is not considered in 
the calculation. 

The calculated results for the unrelax surface energies 
of low index surfaces (100), (110) and (111) are listed in 
Tab.-Ul, in comparison with the relaxed results of Foiles 
et ahEET 



The calculated results for the adsorption energies 
at different sites (see Fig. |^) and the hopping diffu- 
sion barriers t/diff as well as the island formation energies 



on the (100) surface are given in Tab. IV. The result of 
J/diff for Ag is close to the ab initio results 0.52 eV-XLDA) 
and 0.45 eV (GGA) reported by Yu and SchefileiEa. The 
binding energies i^bind of adatom dimers are calculated 
to be negative, suggesting that adatoms tend to form is- 
lands. For trimers on the (100) surface, the calculated 
results disfavor the one dimensional configuration. 
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FIG. 7. Top view of adatoms geometries at (100) sur- 
face, (a) The four-fold hollow site (FFHS); (b) The two-fold 
bridge site (TFBS); (c) One-dimensional (ID) trimer; (d) 
Two-dimensional (2D) trimer. 



TABLE III. The calculated surface energies of low index 
surfaces. The first row are the present results. The second row 
are the theoretical results of Foiles et al.. The experimental 
data (average-face values) as well as the theoretical results 
are all taken from Ref. 27. The units are erg/cm^. 



Surfaces 


Ni 


Pd 


Pt 


Cu 


Ag 


An 


Al 


7100 


1702 


1318 


1485 


1411 


782 


891 


614 




1580 


1370 


1650 


1280 


705 


918 




7110 


1856 


1451 


1650 


1533 


855 


945 


680 




1730 


1490 


1750 


1400 


770 


980 




7111 


1595 


1181 


1286 


1320 


714 


768 


550 




1450 


1220 


1440 


1170 


620 


790 




exp. 


2380 


2000 


2490 


1790 


1240 


1500 





TABLE IV. The calculated adatom adsorption and island 
formation properties on (100) surface. The units are eV. 



Properties 



Ni Pd Pt Cu Ag An Al 



Sad to FFHS -3.59 -3.02 -4.61 -2.77 -2.20 -3.14 -2.85 

Bad to TFBS -2.52 -2.48 -4.10 -2.21 -1.76 -2.93 -2.49 

Uam 1.07 0.54 0.51 0.56 0.44 0.21 0.36 

£;bmd(dimer) -0.42 -0.50 -0.71 -0.37 -0.35 -0.49 -0.29 

_Ebind(lD trimer) -0.82 -0.97 -1.35 -0.71 -0.67 -0.94 -0.55 

_Ebind(2D trimer) -0.90 -1.02 -1.39 -0.78 -0.72 -0.97 -0.58 



D. Molecular dynamics 

In the above subsections the calculations are static. In 
this subsection, the potentials are tested in the constant- 
volume-temperature molecular dynamics (NVT-MD) 
simulation for melting processes for Cu and Pt. The sim- 
ulation box contains 500 atoms. Gear predictor-corrector 
algorithm and Verlet neighbor list are applied. The time 
step is one femtosecond (10~^^ s). The ensemble average 
of the origin-independent translational order parameter 
is calculated after equilibration of at least 5000 steps 



-1 2 1- 

^5^cos(K.R.) -I- l^sin(K.R., 



H 2\ 



(43) 



where A'' is the number of atoms in the box, K is the 
reciprocal basis vector for the initial structure (for ex- 
ample, K = 27r/a(— 1, 1, — 1) for fee lattice), and Ri are 
the position vectors for the atoms. Fig. ^ shows the or- 
der paramters as a function of temperature for Cu and 
Pt. In comparison with experimental data, the melting 
points are underestimated by an amount of 200-300 K_ 
Our result of Cu is worse than that of Foiles and AdamsEHl 
(1340 K), but that of Pt is better than theirs (1480 K). 
Fig.^ shows the pair distribution functions g(r) of Cu at 
different temperatures. 
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FIG. 8. Translational order parameters versus tempera- 
ture for Cu and Pt. The experimental melting points for 
Cu(1358K) and Pt(2045K) are denoted by the two vertical 
solid lines. 
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FIG. 9. Pair distribution functions of Cu at different tem- 
peratures. 

We also simulated the melting of a slab. The size of 
the simulation box is 5x5x20, while that of the slab is 
5x5x10 (containing 21 atomic layers or 1050 atoms). 
The slab is placed at the center of the simulation box, 
which is large enough to ensure that the slab does not 
interact with its images. The atomic configuration is de- 
scribed by the density profile N{z) along the direction 
perpendicular to the slab. N{z) is obtained by averag- 
ing over 1000 steps after running 20000 steps. Surface 
premelting is observed at 900 K, and the liquid fronts 
propagate inward when the temperature rises (920K). At 
950K the slab completely melts. 
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FIG. 10. Density profiles along the direction perpendicular 
to the slab of Cu at different temperatures. 



The simulation results suggest that melting is a 
surface-initiated process. It would be interesting to in- 
vestigate the temperature dependence of the depth of the 
melten layer. On the other hand, does the solidification 
process also begin from the surface, or from the core? 

VI. CONCLUDING REMARKS 

We have presented a systematical method of obtain- 
ing long-range embedded-atom potentials. The model 
parameters can be obtained explicitly from six physical 
inputs and the individual potentials are inverted from the 
analytical functions of lattice sums thus arbitary fitting 
can be avoided. It is shown that the model is able to 
produce satisfactory results of elastic constants, phonon 
eigenfrequencies, phase stabilities, surface properties and 
melting points for the fee transition metals. The poten- 
tials are suitable for computer simulation because of their 
rapid convergence. 

Deriving interatomic potentials from ab initio calcula- 
tions when the experimental data are not available has 
becomej-an abvious trend in the world of material sim- 
ulationHj The reason has been explained very well in a 
recent paper by Payne et al.Eil. In this regard, the present 
method (as well as the method of Bazant and Kaxirasli3) 
may represent an idea of bridging the gap between mate- 
rial theory and electronic structure theory by the method 
of inverting ab initio EAM potentials (or angularly de- 
pendent many-body potentials) from first-principles cal- 
culations. The ab initio binding energy curve can be de- 
composed into repulsive and attractive parts, represent- 
ing the contributions of the pair potential and embedding 
energy respectively. By using our method, the ab initio 
EAM potentials can be obtained by inverting from the 
corresponding parts. 

The present model is easy to be generalized to the al- 
loy case by assuming that the pair potential between 
unlike atoms is given by Johnson's formula VabM = 
{l/2)mr)/fa{r)]VUr) + [fa{r)/ Mr)]V,,ir)}. M Fi- 
nally it should be pointed out that the present model 
fails in predicting the bcc transition metals. Modifying 
the functional forms of p, $ and F{p) does not help much 
to solve this difficulty. 

Acknowledgments: Q.X. gratefully thanks the Max- 
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